﻿using System;

namespace AnthStat.Infrastructure.Statistics
{
    public static class PFromZ
    {
        public static double GetPercentile(double Z)
        {
            // First, convert the Z to
            // type double. 
            double z = Z;
            double y = 0;
            double alnorm = 0;

            // Algorithm AS66 Applied Statistics (1973) vol22 no.3
            // Computes P(Z<z)

            //When z=1.96, the values returned are .975 (upper=false) or .025(upper=true) (roughly).
            //z = parseFloat(z);            
            bool upper = false;

            double ltone = 7.0;
            double utzero = 18.66;
            double con = 1.28;
            double a1 = 0.398942280444;
            double a2 = 0.399903438504;
            double a3 = 5.75885480458;
            double a4 = 29.8213557808;
            double a5 = 2.62433121679;
            double a6 = 48.6959930692;
            double a7 = 5.92885724438;
            double b1 = 0.398942280385;
            double b2 = 3.8052 * Math.Exp(-8.0);
            double b3 = 1.00000615302;
            double b4 = 3.98064794e-4;
            double b5 = 1.986153813664;
            double b6 = 0.151679116635;
            double b7 = 5.29330324926;
            double b8 = 4.8385912808;
            double b9 = 15.1508972451;
            double b10 = 0.742380924027;
            double b11 = 30.789933034;
            double b12 = 3.99019417011;

            if (z < 0)
            {
                upper = !upper;
                z = -z;
            }
            if (z <= ltone || upper && z <= utzero)
            {
                y = 0.5 * z * z;
                if (z > con)
                {
                    alnorm = b1 * Math.Exp(-y) / (z - b2 + b3 / (z + b4 + b5 / (z - b6 + b7 / (z + b8 - b9 / (z + b10 + b11 / (z + b12))))));
                }
                else
                {
                    alnorm = 0.5 - z * (a1 - a2 * y / (y + a3 - a4 / (y + a5 + a6 / (y + a7))));
                }
            }
            else
            {
                alnorm = 0;
            }
            if (!upper) alnorm = 1 - alnorm;

            // After all that math is done, take our decimal and return it to the calling function
            return alnorm * 100;

        } // end GetPercentile
    }
}
